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ABSTRACT 

The masses and radii of low-magnetic field neutron stars can be measured by combining different 
observable quantities obtained from their X-ray spectra during thermonuclear X-ray bursts. One 
I of these quantities is the apparent radius of each neutron star as inferred from the X-ray flux and 

. spectral temperature measured during the cooling tails of bursts, when the thermonuclear flash is 

believed to have engulfed the entire star. In this paper, we analyze 13,095 X-ray spectra of 446 X-ray 
bursts observed from 12 sources in order to assess possible systematic effects in the measurements 
of the apparent radii of neutron stars. We first show that the vast majority of the observed X-ray 
spectra are consistent with blackbody functions to within a few percent. We find that most X- 
ray bursts follow a very well determined correlation between X-ray fiux and temperature, which is 
consistent with the whole neutron-star surface emitting uniformly during the cooling tails. We develop 
rj^ • a Bayesian Gaussian mixture algorithm to measure the apparent radii of the neutron stars in these 

sources, while detecting and excluding a small number of X-ray bursts that show irregular cooling 
behavior. This algorithm also provides us with a quantitative measure of the systematic uncertainties 
' in the measurements. We find that those errors in the spectroscopic determination of neutron-star 

radii that are introduced by systematic effects in the cooling tails of X-ray bursts are in the range 
I , ~ 3 — 8%. Such small errors are adequate to distinguish between different equations of state provided 

O , that uncertainties in the distance to each source and the absolute calibration of X-ray detectors do 

not dominate the error budget. 
^ , Subject headings: stars: neutron — X-rays: bursts 
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CN ! 1. INTRODUCTION 

" The thermal spectra of neutron stars during thermonuclear X-ray bursts have been used during the last three decades 
in numerous attempts to measure the neutron-star masses and radii (e.g., van Paradijs 1978, 1979; Foster, Fabian & 
J — Ross 1986; Sztajno et al. 1987; van Paradijs & Lewin 1987; Damen et al. 1989, 1990). Such efforts were often hampered 
by large systematic uncertainties in the estimated distances to the bursters and in the theoretical models for their 
• . X-ray spectra. Moreover, the relatively small number of X-ray bursts observed by early satellites from each individual 
■ source made it impossible to assess systematic uncertainties related to the degree of anisotropy of the thermonuclear 
burning on the neutron-star surface, or of the obscuration by and the reflection off the accretion flow. 

The situation has changed dramatically in the last few years. The distances to several globular clusters that contain 
bursting neutron stars has been narrowed down with the use of observations by the Hubble space telescope (see, e.g., 

Kuulkers et al. 2003; Ozel, Giiver, & Psaltis 2009; Giiver et al. 2010b). The distances to X-ray bursters in the Galactic 
disk and halo have also been determined using alternate methods (e.g., Giiver et al. 2010a). Theoretical models of 
the X-ray spectra of bursting neutron stars have been greatly improved and can account for the subtle effects of the 
] presence of heavy metals in their atmospheres (e.g., London, Taam, & Howard 1986; Foster, Fabian, & Ross 1986; 
Ebisuzaki 1987; Madej, Joss, & Rozahska 2004; Majczyna et al. 2005; Suleimanov, Poutanen, & Werner 2011). Finally, 
the high signal-to-noise observations of the X-ray spectra of several hundreds of bursts with the Rossi X-ray Timing 
Explorer (Galloway et al. 2008a) allow for the formal uncertainties of individual measurements to be substantially 
reduced. 

The combination of these developments led recently to the measurement of the masses and radii of several neutron 
stars (Ozel et al. 2009; Giiver et al. 2010a, 2010b; Steiner, Lattimer, & Brown 2010) which are already sufficient to 
provide broad constraints on the equation of state of neutron-star matter (Ozel, Baym, & Giiver 2010). In some cases, 
the formal uncertainties in the spectroscopic measurements are as low as a few percent (Giiver et al. 2010b) suggesting 
that their accuracy might be limited by systematic effects rather than by counting statistics. 

In this series of articles, we use the large database of X-ray bursts observed with the RXTE in order to assess 
the systematic effects in various spectroscopic measurements of their properties. In the first article, we focus on the 
measurements of the apparent surface areas of 12 neutron stars as inferred from their X-ray spectra during the cooling 
tails of the bursts. Our goal is to quantify the degree to which: (i) the X-ray burst spectra observed in the RXTE 
energy range can be described by blackbody functions (the so-called color correction arising from atmospheric effects 
are then applied a posteriori); (ii) the entire surface area of each neutron star burns practically uniformly during the 
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References : (1) Harris 1996; (2) Giiver et al. 2010a (3) Juett et al. (2004, 2006); (4) Wroblewski 
et al. 2008; (5) Wijnands et al. 2005; (6) Valenti et al. 2007; (7) Giiver ct al. 2010b; (8) Chevalier 
et al. 1999 

^ Optical/IR observations of the globular cluster 

'^ High resolution spectroscopy of X-ray absorption edges 

'* Average of continuum X-ray spectroscopy 

'^ Optical spectroscopy of the counterpart 

cooling of the btirsts; (Hi) the accretion flows make minor contributions to the emission during the bursts. 

2. SOURCES, OBSERVATIONS, AND DATA ANALYSIS 

2.1. Source and Burst Selection 

We base our study on the X-ray burst catalog of Galloway et al. (2008a). We chose 12 out of the 48 sources in the 
catalog based on the following criteria: 

(i) We considered sources that show at least two bursts with evidence for photospheric radius expansion, based on 
the definition of the latter used by Galloway et al. (2008a). This requirement arises from our ultimate aim, which 
is to measure both the mass and the radius of the neutron star in each system using a combination of spectroscopic 
phenomena (as in, e.g., Ozel et al. 2009). 

(a) We excluded dippers, ADC sources, or known high-inclination sources. This list includes EXO 0748—676, 
MXB 1659-298, 4U 1916-05, GRS 1747-312, 4U 1254-69, and 4U 1710-281, for which it was shown that geo- 
metric effects related to obscuration or reflection significantly affect the flux from the stellar surface that is measured 
by a distant observer (Galloway, Ozel, & Psaltis 2008b). 

(Hi) We did not consider the known mfllisecond pulsars SAX J1808.4-3658 and HETE J1900. 1-2455 because the 
presence of pulsations in their persistent emission implies that their magnetic fields are dynamically important and, 
therefore, may affect the properties of X-ray bursts. 

(iv) We excluded the sources GRS 1741.9-2853 and 2E 1742.9-2929 as weU as a smaU number of bursts from Aql X-1, 
4U 1728—34, and 4U 1746—37, for which there is substantial evidence that their emission is affected by source confusion 
(Galloway et al. 2008a; Keek et al. 2010). 

(v) For each source, we considered only bursts for which the fiux in the persistent emission prior to the burst is at 
most 10% of the inferred Eddington flux for each source, i.e., 7 = i^per/^Edd < 0.1, as calculated by Galloway et al. 
(2008a). Because we subtract the pre-burst persistent emission from the decay spectrum of each X-ray burst, this 
requirement reduces substantially the systematic uncertainties introduced by potential changes in the emission from 
the accretion flow during the burst. 

Table [1] lists all the X-ray bursters that fulfill the above requirements, together with the number of bursts observed 
by RXTE for each source. This is the complete list of sources for which the masses and radii can be measured, in 
principle, using the spectroscopic method of Ozel et al. (2009), with currently available data. Our analyses of the 
bursts from EXO 1745-248, 4U 1608-52, and 4U 1820-30 have been reported elsewhere (Ozel et al. 2009; Guver et 
al. 2010a, 2010b) and will not be repeated here. 

2.2. Data Analysis 

We analyzed the burst data for the 12 sources shown in Table [T] following the method outlined in Galloway et al. 
(2008a; see also Ozel et al. 2009; Giiver et al. 2010a, 2010b). 

For each source, we extracted time resolved 2.5—25.0 keV X-ray spectra using the seextrct ftool for the science event 
mode data and the saextrct ftool for the science array mode data from all of the RXTE/PCA layers. Science mode 
observations provide high count-rate data with a nominal time resolution of 125yLts in 64 spectral channels over the 
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whole energy range (2—60 keV) of the PCA detector. Following Galloway et al. 2008a, we extracted spectra integrated 
over 0.25, 0.5, 1, and 2 s time intervals, depending on the source count rate during the burst, so that the total number 
of counts in each spectrum is roughly constant. (In a few cases, data gaps during the observations result in X-ray 
spectra integrated over shorter exposure times). We took the spectrum over a 16 s time interval prior to the onset of 
each burst as the spectrum of the persistent emission, which we subtracted from the burst spectra as background. 

We generated separate response matrix files for each burst using the PCARSP version 11.7, HEASOFT release 
6.7, and HEASARC's remote calibration database and took into account the offset pointing of the PCA during the 
creation of the response matrix files. This latest version of the PCA response matrix makes the instrument calibration 
self-consistent over the PCA lifetime and yields a normalization of the Crab pulsar that is within 1— tr of the calibration 
measurement of Toor & Seward (1974) for that source. In §4, we discuss in some detail the effect of the uncertainties 
in the absolute flux calibration on the measurement of the apparent surface area of neutron stars. Finally, we corrected 
all of the X-ray spectra for PCA deadtime following the method suggested by the RXTE teani3- 

To analyze the spectra, we used the Interactive Spectral Interpretation System (ISIS), version 1.4.9-55 (Houck & 
Denicola 2000). For each fit, we included a systematic error of 0.5% as suggested by the RXTE cahbration teanfl. 

We fit each spectrum with a blackbody function using the bbodyrad model (as defined in XSPEC, Arnaud 1996) and 
multiplied it by the tbabs model (Wilms, Allen, & McCray 2000) that takes into account the interstellar extinction, 
assuming ISM abundances. The model of the X-ray spectrum for each burst, therefore, depends on only three param- 
eters: the color temperature of the blackbody, Tc, the normalization of the blackbody spectrum. A, and the equivalent 
hydrogen column density, N-^ of the interstellar extinction. 

Allowing the hydrogen column density Nn to be a free parameter in the fitting procedure leads to correlated 
uncertainties between the amount of interstellar extinction and the temperature of the blackbody. Moreover, for 
every 10^^ cm~^ overestimation (underestimation) of the column density, the inferred flux of the thermal emission is 
systematically larger (smaller) by 10~^ erg cm~^ s"""^ (Galloway et al. 2008a). Reducing the uncertainties related 
to the interstellar extinction for each burster is, therefore, very important in controlling systematic effects. 

A recent analysis of high resolution grating spectra from a number of X-ray binaries (Miller, Cackett, & Reis 
2009) showed that the individual photoelectric absorption edges observed in the X-ray spectra do not show significant 
variations with source luminosity or spectral state. This result strongly suggests that the neutral hydrogen column 
density is dominated by absorption in the interstellar medium and does not change on short timescales. In order to 
reduce this systematic uncertainty and given the fact that there is no evidence of variable neutral absorption for each 
system, we fixed the hydrogen column density for each source to a constant value that we obtained in one of the 
following ways, (i) For a number of sources, the equivalent hydrogen column density was inferred independently using 
high-resolution spectrographs, (ii) In cases where only an optical extinction or reddening measurement exists, we used 
the relation given by Giiver & Ozel (2009) to convert it to the equivalent hydrogen column density, (in) Finally, for 
three sources (4U 1702-429, KS 1731-260, SAX J1750.8-2900) there are no independent hydrogen column density 
measurements published in the literature. In these cases, we fit the X-ray spectra obtained during all the X-ray bursts 
of each source allowing the Nh value to vary in a wide range. We then found the resulting mean value of Nh and 
used this as a constant in our second set of fits to all the X-ray bursts. In Table [1] we show the adopted values for 
the hydrogen column density together with the measurement uncertainties, the method with which the values were 
estimated, and the appropriate references. Future observations of these sources with X-ray grating spectrometers 
onboard Chandra and XMM— Newton can help decrease the uncertainty arising from the lack of knowledge of the 
properties of the interstellar matter towards these sources. 

Our goal in this article is to study potential systematic uncertainties in the measurement of the apparent area of 
each neutron star during the cooling tails of X-ray bursts. Hereafter, we adopt the following working definition of 
the cooling tail. It is the time interval during which the inferred flux is lower than the peak flux of the burst or the 
touchdown flux for photospheric radius expansion bursts. For the purpose of this definition, we use as a touchdown 
point the first moment at which the blackbody temperature reaches its highest value and the inferred apparent radius 
is lowest. In order to control the countrate statistics, we also set a lower limit on the thermal flux during each cooling 
tail of 5x10"^ erg s"^ cm~^ (or 5xl0~^° for the exceptionally faint sources 4U 1746-37 and 4U 0513-401). 

In §3 and 4, we discuss in detail our approach of quantifying the systematic uncertainties in the measurements of 
the apparent radii during the cooling tails of thermonuclear bursts, using the sources KS 1731—260, 4U 1728—34, 
and 4U 1636—536 as case studies. For another analysis of the cooling tails of X-ray bursts from 4U 1636—536, see 
Zhang, Mendez, & Altamirano (2011). In §5, we repeat this procedure systematically for seven additional sources from 
Tabled] Finally, the cooling tails of the bursts observed from 4U 0513—40 and Aql X-1 show irregular behavior, and 
we report our analysis of them in the Appendix. 

3. SYSTEMATIC UNCERTAINTIES IN THE SPECTRAL SHAPES 

Our first working hypothesis is that the spectra of neutron stars during the cooling tails of thermonuclear bursts can 
be modeled by blackbody functions in the observed energy range. The fact that X-ray spectra can be described well 
with blackbody functions has been established since the first time resolved X-ray spectral studies of thermonuclear 
bursts (see, e.g., Swank et al. 1977; Lewin, van Paradijs, & Taam 1993; Galloway et al. 2008a and references therein). 
Under that assumption, the temperature measured by fitting blackbodies to the spectra are then corrected for the 

^ ftp://legacy.gsfc.nasa.gov/xte/doc/cook_book/pca_deadtime.ps 

^ http: / / www. universe. nasa.gov / xrays / programs / rxte /pca/doc/ rmf /pcarmf-11.7/ 
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Fig. 1. — The histograms show the distributions of X^/d.o.f. values obtained from fitting 1309, 2519, and 4596 X-ray burst spectra 
observed from the sources KS 1731—260, 4U 1728—34, and 4U 1636—536, respectively. The solid lines show the expected x^/d.o.f. 
distributions for the number of degrees of freedom used during the fits. The vertical dashed lines correspond to the highest values of 
X^/d.o.i. that we considered as statistically acceptable for each source. The vast majority of spectra are well described by blackbody 
functions. 



effects of the atmosphere by applying a numerical prefactor called the color correction factor = T^/Tcs^ which is the 
ratio of the color to the effective temperature. This approach is expected to be only an approximation for a number 
of reasons. 

First, theoretical models of the atmospheres of neutron stars in radiative equilibrium invariably show that the 
emerging spectra are broader than blackbodies, especially at low temperatures and towards low photon energies (e.g., 
Madej et al. 2004; Majczyna et al. 2005). Moreover, for neutron stars that are rapidly spinning, the spectra measured 
by an observer at infinity will be broader and more asymmetric with respect to those at the stellar surface (Ozel & 
Psaltis 2003). The expected degree of broadening and asymmetry will be at most comparable to 

u ^ 2^ ^ / ^3 X / 

c c ■ V600 RzJ \10kmJ ' ^ ' 

where Vg is the spin frequency of the neutron star and R is its radius. This is an uncertainty that can, in principle, be 
corrected for by fitting the X-ray data with theoretical models of the X-ray spectra emerging from the atmospheres of 
rapidly spinning neutron stars. 

Second, the accretion flow in the vicinity of the neutron star may alter in a number of ways the observed spectrum. 
It may act as a mirror, reflecting a fraction of the emerging radiation towards the observer and, therefore, increasing 
the apparent surface area of emission. It might occult a fraction of the stellar surface, reducing the apparent surface 
area of emission (see Galloway et al. 2008b). If the accretion flow is surrounded by a hot corona, Comptonization 
of the surface emission will alter its spectrum as well as the relation between energy flux and apparent surface area 
(e.g., Boutloukos, Miller, & Lamb 2010). Finally, the residual emission of the accretion flow, which we attempt to 
remove by subtracting the pre-burst spectrum of the source, might change during the duration of the burst, introducing 
systematic changes in the net spectrum. 

All of the above effects have the potential of changing the shape of the X-ray spectrum of a burster and, perhaps 
more importantly, lead to stochastic changes in the inferred apparent radii. Their overall effect, however, is expected 
to be significantly reduced by the fact that the intense radiation field during each X-ray burst should either disrupt 
the inner accretion flow or cool any corona to the Compton temperature of the radiation. Such a phenomenon has 
been observed during a supcrburst from 4U 1820—30 (Ballantyne & Strohmayer 2004). Moreover, these effects are 
expected to be the strongest for high inclination systems (as was inferred by, e.g.. Galloway et al. 2008b). This is why 
we not only limited our sample to bursts for which the ratio of the persistent flux prior to the burst to the Eddington 
flux 7 < 10% in order to minimize accretion-related uncertainties but also selected the known high inclination sources 
out of our sample (see §2). 

A final source of systematic uncertainties in the determination of the X-ray spectral shape of a burster is related to 
errors in the response matrix of the detector. These are expected to be < 0.5%, according to the RXTE calibration 
team (see §2.2). 

3.1. Goodness-of-fit Measures 

In order to quantify the degree to which any of these effects influence the spectra of X-ray bursters, we fit all cooling- 
tail spectra from each source with a blackbody model, allowing for a level of systematic uncertainty (Xsys, determined 
in the following way. We reduced the degree of freedom in this procedure by setting dsys in each spectral bin to be a 
constant fraction ^ of the Poisson error (Tformai and fixed the value of ^ to a constant for each source. We then added 
this systematic uncertainty to the formal Poisson error of each measurement in quadrature, i.e., = cr'^o-c-may + '^sys; 
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TABLE 2 







Properties of X'^ 
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0.52 
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40 


b 


2.0 
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4U 1746-37 
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1.6 


85.0% 


SAX J1748.9- 


-2021 


104 


0.26 


1.9 


99.0% 


SAX J1750.8- 


-2900 


82 


0.54 


1.7 


100% 


AQL X-1 




2191 


0.45 


1.6 


92.1% 



^ Maximum value of X^/dof beyond which we consider the blackbody fits of the spectra for each source to be statistically unacceptable. 
^ The distributions for these sources required no addition of systematic uncertainties. 

and defined a statistic X'^ such that 

x'Ks) - YTex' . (2) 

Here, is the standard statistic calculated for each fit. The posterior distribution of the new statistic may not be 
the same as that of the statistic and will depend on the nature of the systematic uncertainties. In order to explore 
the degree and nature of systematic uncertainties, we fit for each source the distribution with the x^ distribution 
expected for the number of degrees of freedom used in the fits, with ^ as a free parameter. Note that the expected 
distribution that we use is formally correct if the countrate data in the individual spectral bins have uncorrelated 
errors. 

The expected x^ a-nd the observed X^ distributions for the optimal ^ value for KS 1731—260, 4U 1728—34, and 
4U 1636—536 are shown in Figure [1] The ^ value required to make the observed X^ distributions for KS 1731—260 
consistent with the expected ones is 0.55 (see Table[2]). This value corresponds to systematic uncertainties that are very 
small. As an example, we consider a typical 0.25 s integration for KS 1731—260, when its spectrum is characterized by 
a color temperature of 2.5 keV. In this case, the average number of counts in each spectral bin in the 3-10 keV energy 
range is 110 ct and the corresponding Poisson uncertainty is ~ 9.5%. Multiplying this by the inferred value ^ — 0.52 
for this source leads to the conclusion that the systematic uncertainties required to render the observed spectrum 
consistent with a blackbody are ~ 5%. For the case of 4U 1728—34, which is brighter than KS 1731—260, the formal 
uncertainties are ~ 6% and the systematic uncertainties amount to ~ 3%. 

Figure [T] shows that the resulting X"^ distributions can be well approximated by the x^ distribution but have weak 
tails extending to higher values. These tails most likely arise from the inconsistency of only a small number of spectra 
with blackbody functions even though we cannot rule out the possibility that the X^ statistic follows a different 
posterior distribution than x^- 

Figure [U allows us also to identify the X-ray spectra that are statistically inconsistent with blackbodies. For each 
source, there is a maximum value of X'^ per degree of freedom beyond which the distribution of X^ values deviates from 
the theoretical expectation. For the case of KS 1731—260, this limiting value of the reduced X^ is 1.7, for 4U 1728—34 
it is 1.7, and for 4U 1636—536 is 1.9 (see Tabled]). The spectra with high reduced X"^ values often occur at the 
late stages of the cooling tails, where the subtraction of the persistent emission is the most problematic. However, 
unacceptable X^ values may also be found in other, seemingly random places of the cooling tails. Nevertheless, the 
fraction of spectra that is inconsistent with blackbodies are < 3%, < 6%, and < 2% for KS 1731—260, for 4U 1728—34, 
and for 4U 1636—536, respectively. Hereafter, we consider only the spectra that we regard to be statistically acceptable, 
given the values of the reduced X^ . 

4. SYSTEMATIC UNCERTAINTIES IN THE INFERRED EMITTING AREA 

Our second working hypothesis is that, during the cooling tail of each burst, the entire neutron star is emitting 
uniformly with negligible lateral temperature variations. This assumption is again expected to be violated at some 
level for a number of reasons. The non-uniformity of accretion onto the neutron star (e.g., Inogamov & Sunyaev 1999), 
the finite time of propagation of the burning front around the star (see, e.g., Nozakura, Ikeuchi, & Fujimoto 1984; 
Bildsten 1995; Spitkovsky, Levin, & Ushomirsky 2002), as well as the excitation of non-radial modes on the stellar 
surface (Heyl 2004; Piro & Bildsten 2005; Narayan & Cooper 2007) are all expected to lead to some variations in the 
effective temperature of emission at different latitudes and longitudes on the stellar surface. 

The characteristics of burst oscillations observed during the cooling tails of X-ray bursts, however, imply that any 
variations in the surface temperatures of neutron stars can only be marginal. Indeed, any component of the variation 
in the surface temperature that is not symmetric with respect to the rotation axis leads to oscillations of the observed 
flux at the spin frequency of the neutron star. Such oscillations have been observed in the tails of bursts from many 
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Fig. 2. — The dependence of the color correction factor fc = T^ft/Tc on the color temperature of the X-ray spectrum for two sequences 
of models by Madej et al. (2004), Majczyna et al. (2005), and Suleimanov et al. (2011) with different metal abundances and surface 
gravities. Between 1 keV and 2.5 keV, the color correction factor depends weakly on the color temperature. 

sources (Galloway et al. 2008a). The r.m.s. amplitudes of burst oscillations in the tails of bursts can be as large as 
15%, although the typical amplitude is ~ 5% (Muno, Ozel, & Chakrabarty 2002). The stringent upper limits on the 
observed amplitudes at the harmonics of the spin frequencies can be accounted for only if the temperature anisotropics 
are dominated by the m = 1 mode in which exactly half the neutron star is hotter than the other half (Muno et al. 

2002) . Moreover, the rather weak dependence of the r.m.s. amplitudes on photon energy (Muno, Ozel, & Chakrabarty 

2003) requires that any temperature variation between the hotter and cooler regions of the neutron star is < 0.2 keV, 
even for the bursts that show the largest burst oscillation amplitudes. All of the above strongly suggest that the 
expected flux anisotropy during the cooling tail of an X-ray burst is < 5 — 10% and, therefore, the expected systematic 
uncertainties in the inferred apparent stellar radius will be half of that value, since the latter scales as the square root 
of the flux. 

A final inherent systematic uncertainty in the spectroscopic measurement of the apparent surface area of a neutron 
star arises from the dependence of the color correction factor on the effective temperature of the atmosphere. In 
Figure [2] we show the predicted evolution of the color correction factor as a function of the color temperature of the 
spectrum as measured by an observer at infinity, for calculations by different groups for different neutron star surface 
gravities and atmospheric compositions (Madej et al. 2004; Majczyna et al. 2005; and Suleimanov et al. 2011). To 
make the models comparable to the observed evolution of the blackbody normalizations presented in the next section, 
we plot the color correction factor against the most directly observed quantity, i.e., the color temperature at infinity. 
When the color temperature at infinity is less than 2.5 keV, purely helium or low metallicity models show — 8% 
evolution of the color correction factor per keV of color temperature at infinity, while above 2.5 keV, they show an 
evolution of 12—20% per keV. In contrast, solar metallicity models show a steady increase with color temperature in 
the 1.5 — 3 keV rangeQ 

We infer the apparent surface area 47ri?^pp of a neutron star by measuring the X-ray fiux i^cooi and the color 
temperature Tc at different intervals during the burst, so that 



aSB(Te//c)4 



(3) 



Here, asB is the Stefan-Boltzmann constant, D is the distance to the source, and fc = Tc/T^g is the color correction 
factor. The value of the color correction factor is dictated by the predominant source of absorption and emission in the 
neutron-star atmosphere. It, therefore, depends on the effective temperature, which determines primarily the ionization 
levels, and on the effective gravitational acceleration, which determines the density profiles of the atmospheric layers. 
Both these quantities decrease as the cooling fiux decreases. If we were to assume a constant value for the color 
correction factor, as it is customary, we would obtain a systematic change in the inferred apparent surface area as 
the cooling flux of the burst decreases with time. Such variations have been discussed in Damen et al. (1989) and in 
Bhattacharyya, Miller, & Galloway (2010). Note that this systematic effect can be corrected, in principle, if the data 
are fit directly with detailed models of neutron-star atmospheres (see, e.g., Majczyna & Madej 2005). 

A potentially important source of uncertainty in the measurements is introduced by the errors in the absolute fiux 
calibration of the RXTE/PCA detector. The current calibration of the PC A and the cross-calibration between X-ray 

^ The color correction factor shows a turnover at different Eddington ratios and correspondingly at different color temperatures depending 
on composition, surface gravity, and gravitational redshift. The rapid evolution presented in Suleimanov et al. (2011) preferentially occurs 
at low Eddington ratios and color temperatures smaller than 1 keV (which we cannot explore observationally in the XTE burst data) and 
for very small surface gravities (logg = 14.0) which correspond only to neutron stars with radii > 15 km. 
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Fig. 3. — (Left) The flux-temperature diagram for aU the spectra in the coohng tails of bursts from (top) KS 1731—260, (Mid- 
dle) 4U 1728—34, and (bottom) 4U 1636—536 that have statisticaUy acceptable values of x'^/d-o.f. The diagonal lines correspond to 
the best-fit blackbody normalization and its uncertainty, as reported in the rightmost column of Table 3. (Right) The distribution of 
measured normalization values of the blackbody spectra in three of the flux intervals we chose. The normalization values for the vast 
majority of spectra fall within a narrowly peaked distribution, with only a number of outliers towards lower (for 4U 1728—34) or higher 
values (for 4U 1636—536) of the normalization. This justifies the assumption that the entire neutron star surface is emitting during the 
cooling tail of a burst witih marginal temperature variations in latitude or longitude. 



satellites have been carried out using the Crab nebula as a standard candle (Jahoda et al. 2006; see also Toor & Seward 
1974; Kirsch et al. 2005; Weisskopf et al. 2010). The uncertainties in the flux calibration can be due to a potential 
overall offset between the inferred and the true flux of the Crab nebula, which may be as large as 10% (Kirsch et al. 
2005; Weisskopf et al. 2010). This can only change the mean value of the inferred apparent area in each source and 
does not alter the observed dispersion. We will explore this issue as well as uncertainties related to the variability of 
the Crab nebula itself (Wilson-Hodge et al. 2011) in more detail in Paper III of this series. 



4.1. Flux- Temperature Diagrams 
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Figure ini (left panels) show the dependence of the emerging flux on color temperature for all the spectra in the cooling 
tails of KS 1731—260, 4U 1728—34, and 4U 1636—536 that we consider to be statistically acceptable (see §3). We chose 
these three sources to use as detailed examples because of the high number spectra obtained for each and the fact that 
they span the range of behavior in cooling tracks that we will discuss below. If the whole neutron star is emitting as 
a single-temperature blackbody and the color correction factor is independent of color temperature, then i^cooi should 
scale as T^. Our aim here is to investigate the systematic uncertainties on the measurement of the apparent surface 
area in each source at different flux levels. We, therefore, divided the data into a number of flux bins and plotted 
in the same figure (right panels) the distribution of the blackbody normalization values for some representative bins. 
The blackbody normalization for each spectrum is defined formally as A = fcooi/csB^<f) although in practice this is 
one of the two measured parameters and the flux is derived from the above definition. According to equation the 
blackbody normalization is equal to A = f~'^R^pp/D^ and we report it in units of (kni/10 kpc)^. If we do not correct 
for the dependence of the color correction factor on the flux, we expect the normalization A to show a dependence on 
color temperature that is the mirror symmetric of that shown in Figure 2. 

The flux-temperature diagrams of KS 1731—260, 4U 1728—34, and 4U 1636—536 share a number of similarities but 
are also distinguished by a number of differences. In all three cases, the vast majority of data points lie along a very 
well defined correlation. This reproducibility of the cooling curves of tens of X-ray bursts per source, combined with 
the lack of large amplitude flux oscillations during cooling tails of bursts, provides the strongest argument that the 
thermal emission engulfs the entire neutron-star surface with very small temperature variations at different latitudes 
and longitudes. 

In 4U 1728—34, a deviation of the data points from the Fcooi ~ correlation is evident at high fluxes and color 
temperatures (Tc > 2.5 keV), which may be due to the evolution of the color correction factor at high Eddington fluxes 
(see the discussion in §4 and Fig. [2]). The same deviation is not evident in KS 1731—260 or 4U 1636—536, for which 
the highest temperatures encountered in the cooling tails were ^2.5 keV. 

Finally, in all three sources, a number of outliers exist at the lowest flux levels, with normalizations that deviate 
from the above correlation. In KS 1731—260 and 4U 1636—536, the outliers correspond to higher normalization values, 
whereas in 4U 1728—34 they correspond to lower normalization values with respect to the majority of the data points. 

Any combination of the effects discussed earlier in this and in the previous section may be responsible for the outliers. 
Non-uniform cooling of the neutron star surface will lead to a reduced inferred value for the apparent surface area. 
Reflection of the surface X-ray emission off a geometrically thin accretion disk will cause an increase in the inferred 
value for the apparent surface area. Finally, Comptonization of the surface emission in a corona may have either effect, 
depending on the Compton temperature of the radiation. 

Our goal in this work is not to understand the origin of the outliers, but rather to ensure that their presence does 
not introduce any biases in the measurements of the apparent surface areas inferred from the vast majority of spectra. 
Indeed, including the outliers in a formal fit of the flux-tcmperature correlation will cause a systematic increase of 
the apparent surface area with decreasing flux in KS 1731—260 and in 4U 1636—536 and a systematic decrease in 
4U 1728—34 (c.f. Damen et al. 1990; Bhattacharyya et al. 2010). This issue was largely avoided in our earlier work 
(Ozel et al. 2009; Giiver et al. 2010a, 2010b) by considering only the relatively early time intervals in the cooling tail 
of each burst, which correspond only to the brightest flux bins. In order to go beyond this limitation here, we consider 
all spectral data for each burst and employ a Bayesian Gaussian mixture algorithm, which is a standard procedure 
for outlier detection in robust statistics (e.g., Titterington, Smith, & Makov 1985; McLachlan & Peel 2000; Huber & 
Ronchetti 2009). 

Our working hypothesis that the main peak of normalizations in each flux interval corresponds to the signal and 
the remainder are outliers is shaped by two aspects of the observations. First, at high fluxes, each histogram can be 
described by a single Gaussian with no evidence or room for a second distribution of what we would call outliers. At 
lower fluxes, when the histograms can be decomposed into two Gaussians, the distribution with the highest peak has 
properties that smoothly connect to those of the single Gaussians at higher fluxes. Second, the Gaussians that we call 
our signal always contain the majority of the data points compared to the distribution of what we call the outliers. 
We take these as our criteria for defining our signal. 

The algorithm, which we describe below in detail, allows us also to measure the degree of intrinsic variation in 
the apparent surface area for each source that is consistent with the width of the flux-temperature correlation. Even 
though we compare spectra in relatively narrow flux bins, the different observing modes as well as the different number 
of PCUs used in each observation result in a wide range of count rates and, hence, in a wide range of formal errors. 
This necessitates the use of the Bayesian approach we describe below as opposed to, e.g., performing fits of the 
distributions over blackbody normalization. 



4.2. Quantifying Systematic Uncertainties with a Bayesian Gaussian Mixture Algorithm for Outlier Detection 

Consider first the situation in which there are no systematic uncertainties and the inferred blackbody normalization 
from every spectrum reflects the apparent surface area of the entire neutron star. If this were the case, the intrinsic 
distribution of blackbody normalizations would be a delta function, while the observed distribution would be broader, 
with a width equal to the average formal uncertainties of the measurements. 

In a more realistic situation, there is an intrinsic range of the inferred blackbody normalizations, caused by a 
combination of the various effects discussed earlier. We assume that, in each flux bin, the intrinsic distribution of 
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Fig. 4. — The histogram shows the distribution of measured blaclcbody normalizations that result from fitting all the available spectra 
during the cooling tails of bursts from 4U 1728—34, when the burst flux was in the range (1 — 2) X 10~* erg cm~^ s~^. The red solid line 
shows the Gaussian mixture model that best describes the data. The black dashed curve shows the underlying Gaussian distribution of 
blackbody normalizations that gives rise to the Gaussian mixture model, when the observational uncertainties and the outliers are taken 
into account. The width of the underlying Gaussian distribution reflects the systematic uncertainty in the measurement of the blackbody 
normalization for this flux interval. 



normalization values is the Gaussian 



Pint (v4|^int, flint) = 



: exp 



lilt 



(4) 



where A^^t is its mean and CTint is its standard deviation. The observed distribution of blackbody normalizations will 
be again broader than the intrinsic distribution because of the formal uncertainties in each measurement. Moreover, 
at low flux levels, a small number of outliers exists, which skews and introduces tails to the observed distribution. 
Our goal here is to quantify the most probable value of the blackbody normalization, i.e., the parameter Aint in the 
distribution as well as the intrinsic range of values in each flux bin, i.e., the parameter CTint, while excluding the 
outliers. 

When the inspection of an observed distribution requires an outlier detection scheme (as, e.g., in the lower right 
panels of Fig. [3]), we model the distribution of outliers by a similar Gaussian, i.e., 



fout(^|^out, CTout) 



1 



V2^ 



: exp 



[A- A, 



2 ) 
out / 



(5) 



with a mean Aout and a standard deviation CTout- Our model distribution function is, therefore, in general the Gaussian 
mixture 

Pmodol(^|^int, CTint, Aut, CTout, V) = -Rnt (^l^int , ^int ) + ?7^out (^| Aut , Cout ) , (6) 

where 77 measures the relative fraction of outliers in our sample. 

We are looking for the parameters in these distributions that maximize the posterior probability P( Aint , CTint , ^out , Cout | data) , 
which measures the likelihood that a particular set of parameters is consistent with the data. Using Bayes' theorem, 
this probability distribution is 



P(Aint, CTint, Aout, CTout, Tyjdata) = CP(data| Ant , CTint , ^out , CTout , ?7)P( Ant )P(CTint )P(^out )P(CTout)P(?7) 



(7) 



where P(data|Aint, CTint, ^out, CTout, v) is the posterior probability that the particular set of normalization data have been 
measured for a given set of parameters and the last five terms are the prior probabilities for the model parameters; 
the constant C ensures that the probability density is normalized. We consider a flat prior probability distribution for 
the central value of each Gaussian in an interval that extends beyond the observed range of normalization values for 
each source. We also consider a flat prior probability distribution for the standard deviation for each Gaussian from 
10~^ (km/10 kpc)^ (which is practically zero) to a value equal to the width of the range of normalizations. Finally, 
we take a flat prior distribution of t] between zero and one. Our results are extremely weakly dependent on the ranges 
of parameter values we considered. 

Given a set of N observed normalization values Ai and their uncertainties at , we calculate the posterior probability 
of the data given a set of model parameters as 



N 

P(data|Aint, CTint, AutO-Qut,??) =^"2]^ dA 



1 



: exp 



(A^Af) 



2a 



[Pint(A|Aint, CTint) +?7Pout(^|Aut, CTout)] • (8) 



Here, Ai is the value of each measurement and ct^ is its corresponding uncertainty, which we assume to be normally 
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Fig. 5. — The dependence of the parameters of the intrinsic distribution of blackbody normalizations on X-ray fiux during the cooling 
tails of thermonuclear X-ray bursts in KS 1731—260, 4U 1728—34, and 4U 1636—536, when the outliers have been removed. Each dot 
represents the most likely centroid of the intrinsic distribution, while each error bar represents its most likely width, as calculated using the 
procedure outlined in §4.2 and depicted in Figure |4] In each panel, the solid and dashed horizontal lines show the best-fit normalization 
and its systematic uncertainty inferred using the flux bins that do not correspond to near-Eddington fluxes and are denoted by flUed circles 
on the error bars. 

distributed; C2 is a normalization constant. Using this posterior probability distribution, we identify the most probable 
values of the model parameters, as well as their uncertainties, following standard procedures. 

Figure 2] shows, as an example, the application of the Gaussian mixture algorithm to the distribution of blackbody 
normalizations obtained from fitting the spectra of 4U 1728—34, when the burst flux was in the range (1 — 2) x 
10~* erg cm~^ s~^. The histogram shows the distribution of blackbody normalizations measured in this flux inter- 
val. The parameters of the Gaussian mixture that maximize the posterior probability distribution calculated using 
equation (O correspond to an intrinsic distribution with a mean and standard deviation of Ai^t = 126.1 (km/10 kpc)^ 
and CTint — 13.4 (km/10 kpc)^, respectively, and a distribution of outliers with a mean of Aout = 62.3 (km/10 kpc)^, 
a standard deviation of Uout — 6.1 (km/10 kpc)^, and a relative normalization of 77 = 0.2. In order to compare the 
Gaussian mixture with the observed histogram, we first convolve it with a third Gaussian distribution with a standard 
deviation of (Tformai = 16.5 (km/10 kpc)^ that is equal to the average formal errors of all measurements in this fiux 
interval. The result is shown as a red solid line in Figure SI which demonstrates that our assumption of a Gaussian 
distribution for the majority of the data and for the outliers is consistent with the observations. 

We can appreciate the importance of excluding the outliers from the sample using the Bayesian Gaussian mixture 
algorithm in the following way. If we repeat the above procedure for the same fiux interval in 4U 1728—34 but do not 
allow for the presence of outliers, the parameters of the intrinsic distribution will be Ant = 118.6 (km/10 kpc)^ and 
flout — 26.4 (km/10 kpc)^, which are significantly different than the results quoted above. Moreover, had we assumed 
that there are no systematic errors but that the underlying distribution was a delta function, the most probable 
value would have been ^int = 117.llgg (km/10 kpc)^. The formal errors in such a measurement would have been 
substantially smaller compared to the systematic uncertainties. 

In the following, we show in detail the application of this Gaussian mixture algorithm to all the fiux intervals for 
the three sources. 



4.3. Dependence of Blackbody Normalizations on X-ray Flux 

Figure [5] shows the result of the outlier detection procedure as applied to the various flux intervals of KS 1731—260, 
4U 1728—35, and 4U 1636—536 (see also Table Each dot represents the most likely centroid of the intrinsic 
distribution of blackbody normalizations, while each error bar represents its most likely width. 

The normalization of the blackbody in the case of 4U 1728—34 shows a strong dependence on flux when the source is 
emitting at near-Eddington rates. This is not seen in KS 1731—260 and 4U 1636—536, for which the color temperature 
never reached temperatures above 2.5 keV. In fact, in our observed sample, the strong dependence of the normalization 
on the flux is seen in all three sources with color temperatures in excess of 2.5 keV, namely 4U 1728—34, 4U 1702—429, 
and 4U 1724—307 (see Fig. 6 below) but is absent from the other sources. The evolution of the blackbody normalization 
at near-Eddington fluxes depends on the local gravitational acceleration {g ^ M/R^) on the neutron-star surface and 
on its composition. Therefore, modeling the decline of the blackbody normalization at large fluxes would allow us, in 
principle, to measure the combination Mj of the mass and radius of the neutron star, as has been attempted by 
Majczyna & Madej (2005) and by Suleimanov et al. (2011). 

At low burst fluxes, the observed normalizations show a small albeit statistically signiflcant trend towards lower 
values. In the case of 4U 1728—34, where the potential decrease is the largest, the normalization changed from 
133.7 ± 16.4 (km/lO kpc)^ to 114.0 ± 15.6 (km/10 kpc)^ as the flux declined from 6 x 10~*^ erg s'^ cm'^ to 0.5 x 
10~^ erg s~^ cm~^. This corresponds to a < 15% reduction in the normalization. In other sources, the decline is 
even smaller. This weak dependence seen in the data argues against the models with solar composition, as shown in 
Figure [3 

The decline in the normalization can, in principle, be accounted for with a < 4% increase in the color correction 
factor (see eq. [3]) towards low temperatures. As can be seen in Figure [2l current atmosphere models do not predict 
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TABLE 3 

Blackbody Normalizations'' 













Flux interval (10 


erg s 


cm ^) 










Source 




0.5-1 


1-2 


2-3 


3-4 


4-5 


5-6 


6-7 


7-8 


8-9 


9-10 


Average'' 


4U 1636-536 




112.6 


124.6 


126.9 


142.7 


144.2 


145.7 


145.7 


149.5 






130.7 






±18.7 


±16.1 


±12.4 


±22.4 


±18.4 


±15.1 


±14.1 


±4.6 






±20.9 


4U 1702-429 




151.0 


167.6 


180.4 


184.2 


185.7 


171.4 


151.0 


119.3 


98.2 




176.6 






±13.9 


±10.1 


±9.9 


±4.1 


±3.9 


±8.1 


±8.1 


±17.2 


±12.4 




±11.6 


4U 1705-44 




80.2 
±9.9 


86.9 
±7.1 


86.9 
±10.9 


82.4 
±7.4 














83.9 
±9.1 


4U 1724-307 




98.7 


108.3 


120.4 


127.4 


125.4 


102.8 










113.8 






±7.1 


±4.6 


±12.6 


±3.7 


±19.4 


±15.4 










±15.4 


4U 1728-34 




114.1 


126.1 


137.4 


145.0 


135.2 


137.4 


133.7 


126.1 


117.8 


93.7 


134.4 






±15.6 


±13.4 


±10.9 


±8.9 


±18.2 


±18.4 


±16.4 


±13.1 


±15.4 


±13.1 


±14.9 


KS 1731-260 




72.6 
±6.1 


89.9 
±5.4 


93.0 
±4.1 


95.2 
±9.6 


95.2 
±1.9 












96.0 
±7.9 


4U 1735-44 




71.6 
±4.6 


71.6^= 

+2.0 
-1.5 


72.6<= 

+2.4 
-1.9 
















72.1': 

+ 1.3 
-1.0 


4U 1746-37"' 




12.9 


16.3 


17.3 


15.6 


15.2 


15.2 


19.9 








15.7 






±1.9 


±0.5 


±0.4 


±2.4 


±2.4 


±1.6 


±6.4 








±2.4 


SAX J1748.9- 


-2021 


92.7 
±11.9 


87.7 
±6.6 


91.7 
±7.6 


87.7 
±15.4 














89.7 
±9.6 


SAX J1750.8- 


-2900 


126.9 
±40.5 


110.8 
±19.9 


106.8 
±10.1 


97.2 
±8.1 


86.2 
±8.6 












93.2': 
±9.4 



The parameters of the intrinsic distribution of blackbody normalizations in different flux intervals for all the sources we considered in 
this manuscript. All normalizations are in units of (km/10 kpe)^. 
^ Taking into account all flux intervals for which the color temperature of the spectrum is < 2.5 keV; see text for details. 

The range of normalizations in this flux interval are consistent with no systematic uncertainties; the quoted uncertainties are purely 
statistical. 

For 4U 1746—37, all flux intervals are in units of 10~^ erg s~' cm~^. 



TABLE 4 

The Apparent Radii of X-ray Bursters 



Source 






Rfcm / Diofep, 


4U 1636- 


-53 




11.4 ± 1.0 


4U 1702- 


-429 




13.3 ± 0.4 


4U 1705- 


-44 




9.2 ± 0.5 


4U 1724- 


-307 




10.7 ± 0.7 


4U 1728- 


-34 




11.6 ± 0.7 


KS 1731- 


-260 




9.8 ± 0.4 


4U 1735- 


-44 




o c+O.OS 
°-°-0.06 


4U 1746- 


-37 




4.0 ± 0.3 


SAX J1748.9- 


-2021 


9.5 ± 0.5 


SAX J1750.8- 


-2900 


9.6 ± 0.5 



such an evolution. However, as we discussed in the beginning of §4, a number of effects related to the physics of 
bursts on the neutron star surface (uneven cooling and burst oscillations), as well as the reduced sensitivity of the 
PCA at low energies, are capable of causing the observed trend in the normalization. Moreover, the < 4% effect in 
the data we would aim to model is comparable to the uncertainty in inferring theoretically the color correction factor 
from the atmosphere models and is also comparable to the deviation of the observed and theoretical spectra from 
blackbodies. The latter concern can be remedied by fitting directly theoretical model spectra to the data, but reducing 
the theoretical uncertainties present in the models to less than few percent is significantly more challenging. 

An alternative approach is to allow for a range of values for the color correction factor that span the spread of 
theoretical uncertainties, metallicities, and fluxes of the cooling tails. In earlier work (Giiver et al. 2010a, 2010b) we 
allowed for a 7% range in the color correction factor (from 1.3 to 1.4), which is adequate to account for the evolution in 
the blackbody normalizations shown by the data. Following this approach, we need to identify the range of blackbody 
normalizations we obtained from the data at different flux intervals as a systematic uncertainty in the measurement. 
In order to achieve this, we applied the Bayesian Gaussian mixture algorithm to the combined data set for each source 
for all flux intervals that correspond to a color temperature < 2.5 keV. This way, we computed the most probable 
value for the blackbody normalization for each source in a wide flux range as well as the systematic uncertainties in 
that measurement, which account for the potential decline of the normalization with decreasing flux. 

In Figure [5l we identify the flux intervals we used for each source with a filled circle in the middle of each error 
bar. We also depict with a solid line the most probable value of the normalization in this wide flux range and 
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Fig. 6. — (Left) The distribution of the X^/dof values obtained from fitting the spectra during the tails of thermonuclear X-ray bursts 
observed from the sources 4U 1702—429, 4U 1705—44, and 4U 1724—307 together with the theoretically expected distribution; the vertical 
dashed line marks the maximum value of Jf^/dof beyond which we consider the blackbody model to be inconsistent with the data. 
(Middle) The flux-temperature diagrams of the cooling tails of bursts from the same sources; the solid and dashed line correspond to the 
most probable values of the blackbody normalizations throughout the bursts and their systematic uncertainties. (Right) The dependence 
of the parameters of the intrinsic blackbody normalization on X-ray flux; the solid and dashed lines correspond to the most probable values 
of the normalizations and their systematic uncertainties for the flux bins that are marked by a filled circle. Error bars without filled circles 
appear at near-Eddington fluxes where the color correction factor increases, causing the apparent decline in the normalization. 

with dashed hnes the range of systematic uncertainties. Our results are summarized in Table 3. The blackbody 
normalizations for KS 1731-260, 4U 1728-34, and 4U 1636-536 are 96.0±7.9 (km/10 kpc)^, 134.4±14.9 (km/10 kpc)^ 
and 130.7±20.9 (km/10 kpc)^, respectively, with the uncertainties dominated entirely by systematics. 



5. THE APPARENT RADII OF X-RAY BURSTERS 

Figures [MS] and Tables [2H3] show the results obtained after applying the procedure outlined in §4 to seven more 
sources from Table [T] The last two sources, 4U 0513—40 and Aql X-1 show large variations in the cooling tails of 
individual bursts and are discussed in detail in the Appendix. Moreover, in the Appendix we also discuss a number of 
bursts from 4U 1702—429 and a long burst from 4U 1724—307, which we did not include in the analysis. 

As in the case of the three sources discussed in the previous section, the vast majority of the spectra observed during 
the cooling tails of X-ray bursts are very well described by blackbody functions, with only marginal allowance for 
systematic variations. In the flux-temperature diagrams, the cooling tails largely follow a well defined and reproducible 
track. Finally, in each flux interval of most sources, the range of blackbody normalizations obtained from a large number 
of bursts is consistent with a small degree of systematic uncertainties. 

The measurements of the average blackbody normalizations in all sources are dominated by systematic uncertainties 
(the main exception is 4U 1735—44). However, these uncertainties are small, ranging from ~ 6% in the case of 
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Fig. 7.— Same as Figure|6]but for the sources 4U 1735-44, 4U 1746-37 and SAX J1748.9-2021. 
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Fig. 8. — Same as Figure [6l but for the source SAX J1750.8— 2900. Note that, for this source, the low number of spectra precluded the 
clear identification of outliers in the individual flux intervals and resulted in large systematic uncertainties. However, when all flux intervals 
were considered together, a population of outliers (coming primarily from one burst) became clear and the intrinsic distribution was found 
to be consistent with having no systematic uncertainties. 
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Fig. 9. — (Left) The distribution of x^/dof obtained from fitting the spectra of 40 weak X-ray bursts observed from 4U 1702—429; the 
solid line shows the expected distribution for the same number of degrees of freedom. (Right) The relation between the inferred blackbody 
normalization and the x^/dof observed during the cooling tails of the 40 weak bursts observed from 4U 1702—429. 
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Fig. 10. — (Left) The distribution of x^/dof obtained from fitting the spectra in the cooling tail of one long X-ray burst observed from 
4U 1724—307; the solid line shows the expected distribution for the same number of degrees of freedom. (Right) The flux-temperature 
diagram of the long burst observed from 4U 1724—307 (filled circles) compared to that of the short bursts discussed in the main text (open 
circles). 

4U 1702—429 to ~ 16% in the case of 4U 1636—536. The apparent radius of each neutron star scales as the square 
root of the blackbody normalization. As a result, the errors in the spectroscopic determination of neutron-star radii 
introduced by systematic effects in the cooling tails of X-ray bursts are in the range ~ 3 — 8% (see Table S]). Such 
small errors by themselves do not preclude distinguishing between different equations of state of neutron-star matter. 
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APPENDIX 
A. 4U 1702-429 

A total of 46 bursts have been observed from the source 4U 1702—429. Six of these bursts reach high fluxes and 
have been categorized as Photospheric Radius Expansion bursts by Galloway et al. (2008a). The remaining 40 bursts 
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typically reach lower fluxes. In the main body of this paper, we focused on the 6 bright bursts from this source and 
exclude the remaining for two reasons that we explore in this appendix. 

The left panel of Figure |9] shows the histogram of x^/dof values for the 880 spectra observed during the cooling tails 
of the 40 bursts and compares them to the expected distribution given the number of degrees of freedom. It is evident 
that blackbody functions provide statistically unacceptable fits to the majority of the spectra. 

Had we not considered these spectra unacceptable, we would find that, approximately 4 seconds after the bursts 
start, the inferred blackbody normalizations increase rapidly to large values. The increase in normalization is, in fact, 
correlated with an increase in the x^/dof values, as shown in the right panel of Figure |9l rendering them even more 
suspect. These two arguments strongly suggest that the spectra of the 40 weak X-ray bursts from 4U 1702—429 are 
not dominated by the thermal emission from the neutron star. 

B. 4U 1724-307 

The first X-ray burst observed from 4U 1724—307 on 1996 November 08 during observation 10090-01-01-02 is very 
different compared to the other two X-ray bursts that we used in the main text. The first burst is unusually long 
(~ 120 s), while the other two bursts are much shorter (~ 15 — 20 s range, as is typical for the vast majority of 
the Type-I bursts observed). Moreover, while the spectra of the short bursts are accurately modeled with blackbody 
functions (see top panels of Fig.[6|), fitting the spectra during the cooling tail of the long burst results in unacceptably 
large values of x^/dof and the distritubiton of the resulting x^/dof values do not follow the expected x^ distribution 
(see left panel of Fig. [TU|) . Adding a certain amount of systematic uncertainty to the data can in principle result in 
a decrease in the individual x^/dof values; however, it can not change the fact that the resulting x^/dof distribution 
does not follow the expected distribution. 

If we went ahead and further analyzed these statistically unacceptable spectra, we would have obtained values for 
the blackbody normalization that are larger compared to those of the shorter bursts (see right panel of Fig. [TO)) . 

The large values and the seemingly random distribution of x^/dof suggest that the X-ray spectra in the cooling tail 
of the long burst from 4U 1724—307 are not dominated by the thermal emission from the neutron star. In fact, in't 
Zand & Weinberg (2010) found evidence for atomic edges and a reflection component in the spectra of this burst, 
which they attributed to the presence of heavy metals in the ashes of previous bursts that were exposed by the long 
radius-expansion episode. For these reasons, we excluded the long X-ray burst observed from 4U 1724—307 from our 
analysis. Note that Suleimanov et al. (2011) chose the spectra from this burst in order to infer the mass and radius of 
the neutron star in 4U 1724—307. Had they chosen the other two bursts from the same source, the spectra of which are 
actually well described by blackbody functions, they would have obtained a radius of the neutron star that is ~ 30% 
smaller, making their result consistent with the radii inferred for the other sources using this method (Ozel et al. 2009; 
Giiver et al. 2010a, 2010b) as well as with radii inferred from quiescent neutron stars in globular clusters (e.g., Webb 
& Barret 2007; Guillot, Rutledge, & Brown 2010). 

C. 4U 0513-40 

Six bursts have been observed from the source 4U 0513—40. While most of the X-ray spectra during the cooling tails 
of these bursts can be described well by blackbody functions (left panel of Fig. [TT|) . the fiux-temperature diagrams 
show irregular behaviour for the cooling tails (right panel of Fig. Ill[) . In particular, two of the bursts show a cooling 
tail that is reminiscent of the long burst from 4U 1724—307 discussed above: early in the cooling phase the temperature 
of the blackbody decreases while the flux remains constant and, when the fiux starts decreasing, the cooling occurs 
with blackbody normalizations that are on average larger compared to the other bursts. 

It is plausible that the spectra of these bursts from 4U 0513—40 are also affected by the presence of atomic lines, 
but the relatively low flux of 4U 0513—40 (which is a factor of ~ 3 — 5 smaller compared to that of 4U 1724—307 at 
similar temperatures) still allows us to fit the observed spectra with blackbody functions. We do not consider this 
source suitable for a radius measurement based on the burst data. 

D. AQL X-l 

The transient source Aql X-l has shown 51 X-ray bursts, from which we extracted 2191 spectra. Most of these 
spectra are well fit by blackbody functions, but their cooling tracks in the flux-temperature diagram depend strongly 
on the properties of the bursts (see Fig. [T^ . In particular, a large number of bursts are relatively short and reach 
modest fluxes (typically < 5 x 10~^ erg s~^ cm~^) with cooling tracks in the flux-temperature diagram that are 
reproducible. On the other hand, many bursts, which are relatively longer, reach fluxes that are a factor of ~ 2 higher, 
are often characterized by radius expansion episodes, and follow a range of cooling tracks in the fiux-temperature 
diagram with blackbody normalizations that are typically larger than those of the other bursts. 

One plausible exlanation is related to the fraction of the neutron star surface that is engulfed by the thermonuclear 
flash. In weak flashes, the burning front may not propage across the entire surface and, therefore, only a fraction of 
the stellar surface contributes to the normalization of the blackbody. On the other hand, in stronger flashes, a larger 
fraction of the neutron star surface is covered. 

The source Aql X-l is one of the two main examples (the other being the source 4U 1636—536) discussed by 
Bhattacharrya et al. (2010) as evidence for systematic variations in the burning area during thermonuclear X-ray 
bursts. It is important to emphasize here, however, that the case of Aql X-l is very unusual compared to the remaining 
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Fig. 11. — Same as Figure|6]but for the source 4U 0513-401. 
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Fig. 12. — Same as Figure|6]but for the source Aql X-1. 

sources that we analyzed and is, perhaps, related to the presence of a weak but dynamically important magnetic field 
in this source, as inferred from the observation of intermittent persistent pulsations (Casella et al. 2008). 
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